#include <iostream>
#include <cstdio>

#include "TopologicalOperation.h"
#include "TopologicalModel.h"
#include "Tetrahedron.h"
#include "Node.h"
#include "Setting.h"
#include "FEMModel.h"

#include "Collapse.h"

/*
Do rivara refinement, split the longest edge of a tet over the midpoint.

MC RIvara and C LEvin a 3d refinement algorithm suitable for
adaptive and multi-grid techniques.

Comm. in Appl. Numer. Meth. 8 (1992) 281-290

*****

Apparently, it doesn't work in "real life." It takes lots of CPU,
though. Or maybe I misunderstood the idea of this technique. 

*/
void SplitLongestEdge (FiniteElementModel*fem, Tetrahedron* t)
{
	cout << "Splitting tetra " << t->Number ()  << " ... " <<flush;
	Simplex max;
	Real maxlength = 0.0;
	for (int i=4; i--;)
		for (int j = i; j--;)
		{
			Node *n[] = {t->GetNode (i), t->GetNode (j)};

			Real l = (n[0]->Location (false) - n[1]->Location (false)).Length ();

			if (l >maxlength)
			{
				max = Simplex (EDGE, n, 0);
				maxlength = l;
			}
		}
		TrySplitEdge (fem, max, t, Centroid (max, false));
		cout << "done"<<endl;
}

void ImproveMesh (FiniteElementModel*fem)
{
	Real max = GetNumberSetting ("collapse-degeneracy-passes");
	cout<< "Improving mesh ..." <<flush; 


	int last_degens_detected = 0;
	int iter = 0;

	Real aspect_ratio_threshold = GetNumberSetting ("improvement-aspect-ratio-threshold");


	Real last_mean_ar=0.0;
	do
	{
		int degens_detected =0;
		if (iter >= max)
			break ;

		iter ++;
		if (last_degens_detected != 0)
			cout << " found " << last_degens_detected << " degens \n"
			<< "******* trying another time ********\n\n";

		LinkArray<Tetrahedron> todo (fem->TetrahedronArray ());


		Real total_as_ratio = 0.;
		for (int i =todo.Size (); i--; )
		{
			Tetrahedron* tet = todo[i];
			if (!tet->Valid ())
				continue;

			Real a  = AspectRatio (tet, false);
			if (a > aspect_ratio_threshold)
				continue;

			total_as_ratio += a;

			cout << "tet "  << tet->Number () << "bad: aspect "<< a << flush;

			degens_detected ++;
			SplitLongestEdge (fem, tet);
		}

		Real mean_ar = total_as_ratio / degens_detected;
		if (last_mean_ar > mean_ar)
		{
			cout << "uhoh. Decreasing mesh quality!"<<endl;
			cout << " now: " << mean_ar << " prev run " << last_mean_ar <<endl;
			return;
		}
		last_mean_ar = mean_ar;


		if (degens_detected == 0)
			break;

		if ( last_degens_detected && degens_detected >= last_degens_detected)
		{
			cout << "uhoh. Increasing # of degens!"<<endl;
			return ;
		}


		last_degens_detected = degens_detected ;
	}
	while (1);



	printf ("done\n");
}
